{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "93021c58",
   "metadata": {},
   "source": [
    "## 代码说明  \n",
    "该代码用于计算NDVI的相应统计值，比如mean, variance, std, cv等"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aaff9ccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.windows import Window\n",
    "from rasterio.plot import show\n",
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "from ipywidgets import Dropdown, interact\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f158c45b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# === 1. 初始化输入 ===\n",
    "input_files = glob.glob('D:/0Datasets/NDVI_MAX_MODIS/*.tif')  # 替换为实际路径\n",
    "output_path = \"D:\\\\0GPPvsWater\\\\Modis\\\\统计量\"\n",
    "\n",
    "# 确保文件按年份排序\n",
    "input_files.sort()\n",
    "\n",
    "# === 2. 分块参数 ===\n",
    "CHUNK_SIZE = 256  # 根据内存调整"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a2e35bca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating global mask...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9f0d32dc082b4f4eae5854daeb97fb22",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/25 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# === 3. 获取元数据及初始化全局掩码 ===\n",
    "with rasterio.open(input_files[0]) as src:\n",
    "    meta = src.meta.copy()\n",
    "    height, width = src.height, src.width\n",
    "    nodata = src.nodata\n",
    "    meta.update(dtype=rasterio.float32, count=1, nodata=nodata)\n",
    "\n",
    "# 初始化全局掩码（所有年份的nodata并集）\n",
    "global_mask = np.zeros((height, width), dtype=bool)\n",
    "\n",
    "# === 4. 生成全局掩码（所有年份的无效区域）===\n",
    "print(\"Generating global mask...\")\n",
    "for file in tqdm(input_files):\n",
    "    with rasterio.open(file) as src:\n",
    "        for i in range(0, height, CHUNK_SIZE):\n",
    "            for j in range(0, width, CHUNK_SIZE):\n",
    "                win = Window(j, i, min(CHUNK_SIZE, width - j), min(CHUNK_SIZE, height - i))\n",
    "                chunk = src.read(1, window=win)\n",
    "                # 合并当前文件的无效区域到全局掩码\n",
    "                global_mask[i:i+win.height, j:j+win.width] |= (chunk == nodata) | np.isnan(chunk)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f721690c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating statistics...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "38a3ff6a5ae044e7870b7fbbe5dd80cb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/25 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# === 5. 计算统计量（仅处理有效区域）===\n",
    "sum_ndvi = np.zeros((height, width), dtype=np.float64)\n",
    "sum_sq_ndvi = np.zeros((height, width), dtype=np.float64)\n",
    "valid_counts = np.zeros((height, width), dtype=np.uint16)  # 有效值的个数，用于计算均值等\n",
    "\n",
    "print(\"Calculating statistics...\")\n",
    "for file in tqdm(input_files):\n",
    "    with rasterio.open(file) as src:\n",
    "        for i in range(0, height, CHUNK_SIZE):\n",
    "            for j in range(0, width, CHUNK_SIZE):\n",
    "                win = Window(j, i, min(CHUNK_SIZE, width - j), min(CHUNK_SIZE, height - i))\n",
    "                chunk = src.read(1, window=win)\n",
    "                # 仅处理非全局掩码区域\n",
    "                local_mask = global_mask[i:i+win.height, j:j+win.width]\n",
    "                valid_pixels = (chunk != nodata) & (~np.isnan(chunk)) & (~local_mask)\n",
    "                \n",
    "                sum_ndvi[i:i+win.height, j:j+win.width] += np.where(valid_pixels, chunk, 0)\n",
    "                sum_sq_ndvi[i:i+win.height, j:j+win.width] += np.where(valid_pixels, chunk**2, 0)  # 计算方差，使用下方的公式\n",
    "                #　方差 = E[X²] - (E[X])² = (sum_sq_ndvi / valid_counts) - (sum_ndvi / valid_counts)²\n",
    "                valid_counts[i:i+win.height, j:j+win.width] += valid_pixels.astype(np.uint16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fbe20395",
   "metadata": {},
   "outputs": [],
   "source": [
    "# === 6. 计算结果（修复sqrt错误）===\n",
    "def masked_divide(a, b):\n",
    "    result = np.full_like(a, nodata, dtype=np.float32)\n",
    "    valid = (b > 0) & (~global_mask)\n",
    "    np.divide(a, b, out=result, where=valid)\n",
    "    return result\n",
    "\n",
    "mean_ndvi = masked_divide(sum_ndvi, valid_counts)\n",
    "# 计算方差，　方差 = E[X²] - (E[X])² = (sum_sq_ndvi / valid_counts) - (sum_ndvi / valid_counts)²\n",
    "variance_ndvi = masked_divide(sum_sq_ndvi, valid_counts) - mean_ndvi**2\n",
    "\n",
    "# 修复关键错误：对variance_ndvi中的负值进行保护\n",
    "variance_ndvi = np.where(variance_ndvi < 0, 0, variance_ndvi)  # 方差不应为负\n",
    "std_ndvi = np.full_like(variance_ndvi, nodata, dtype=np.float32)\n",
    "np.sqrt(variance_ndvi, out=std_ndvi, where=(~global_mask))  # 仅对有效区域计算平方根\n",
    "\n",
    "cv_ndvi = masked_divide(std_ndvi, mean_ndvi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bda241e3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved: D:\\0GPPvsWater\\Modis\\统计量/ndvi_mean_2000_2024.tif\n",
      "Saved: D:\\0GPPvsWater\\Modis\\统计量/ndvi_variance_2000_2024.tif\n",
      "Saved: D:\\0GPPvsWater\\Modis\\统计量/ndvi_std_2000_2024.tif\n",
      "Saved: D:\\0GPPvsWater\\Modis\\统计量/ndvi_cv_2000_2024.tif\n"
     ]
    }
   ],
   "source": [
    "# === 7. 保存结果 ===\n",
    "output_stats = {\n",
    "    'mean': mean_ndvi,\n",
    "    'variance': variance_ndvi,\n",
    "    'std': std_ndvi,\n",
    "    'cv': cv_ndvi\n",
    "}\n",
    "\n",
    "for stat_name, stat_data in output_stats.items():\n",
    "    output_file = output_path + f'/ndvi_{stat_name}_2000_2024.tif'\n",
    "    with rasterio.open(output_file, 'w', **meta) as dst:\n",
    "        dst.write(stat_data, 1)\n",
    "    print(f\"Saved: {output_file}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9107f52e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cfb2d5f005134a748d8bb95e6f0abc58",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Statistic:', options=('mean', 'variance', 'std', 'cv'), value='mea…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_stat(stat_name)>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# # === 8. Jupyter交互显示 ===\n",
    "# def plot_stat(stat_name, output_path):\n",
    "#     output_file = output_path + \"/\"  + f'ndvi_{stat_name}_2000_2024.tif'\n",
    "#     plt.figure(figsize=(12, 8))\n",
    "#     with rasterio.open(output_file) as src:\n",
    "#         show(src, title=f'NDVI {stat_name} (2000-2024)')\n",
    "#         plt.colorbar(label='Value')\n",
    "#     plt.show()\n",
    "\n",
    "# interact(plot_stat, stat_name=Dropdown(options=list(output_stats.keys()), description='Statistic:'),output_path=output_path)\n",
    "# === 8. Jupyter交互显示 ===\n",
    "import os\n",
    "\n",
    "def plot_stat(stat_name):\n",
    "    # 在函数内部直接使用外部变量 output_path\n",
    "    output_file = os.path.join(output_path, f'ndvi_{stat_name}_2000_2024.tif')\n",
    "    \n",
    "    if not os.path.exists(output_file):\n",
    "        print(f\"文件不存在: {output_file}\")\n",
    "        return\n",
    "    \n",
    "    # 创建图形和坐标轴\n",
    "    fig, ax = plt.subplots(figsize=(12, 8))\n",
    "    \n",
    "    with rasterio.open(output_file) as src:\n",
    "        # 读取数据\n",
    "        data = src.read(1)\n",
    "        \n",
    "        # 创建掩膜，排除nodata值\n",
    "        if src.nodata is not None:\n",
    "            masked_data = np.ma.masked_where(data == src.nodata, data)\n",
    "        else:\n",
    "            masked_data = data\n",
    "        \n",
    "        # 显示图像并保存返回的图像对象\n",
    "        im = ax.imshow(masked_data, cmap='viridis')\n",
    "        ax.set_title(f'NDVI {stat_name} (2000-2024)')\n",
    "        \n",
    "        # 添加colorbar\n",
    "        plt.colorbar(im, ax=ax, label='Value')\n",
    "        \n",
    "        # 移除坐标轴刻度\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "# 创建交互式控件（移除了output_path参数）\n",
    "interact(\n",
    "    plot_stat, \n",
    "    stat_name=Dropdown(\n",
    "        options=list(output_stats.keys()), \n",
    "        description='Statistic:',\n",
    "        value=list(output_stats.keys())[0]  # 设置默认值\n",
    "    )\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
